Demonstration: Compare Density Models

This notebook demonstrates how to compare the GEODYN output for different density models

First, load the data using the reading functions.

[1]:
%load_ext autoreload
%autoreload 2

import pandas as pd
import numpy as np

##################################
# INPUT PARAMETERS:
##################################
sat = 'st'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001
Accel_Status = 'acceloff'

##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/'+sat+'/results/'+msis86_model+'/'+  msis86_model+'_'+ Accel_Status +'/'

dtm87_model = 'dtm87'
path_to_dtm87 = '/data/runs_geodyn/'+sat+'/results/'+dtm87_model+'/'+  dtm87_model+'_'+ Accel_Status +'/'

jaachia71_model = 'jaachia71'
path_to_jaachia71 = '/data/runs_geodyn/'+sat+'/results/'+jaachia71_model+'/'+  jaachia71_model+'_'+ Accel_Status +'/'


[2]:
import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_read_geodyn_output/')
from a_ReadGEODYN import Read_GEODYN_func

(AdjustedParams_msis86,
 Trajectory_msis86,
 Density_msis86,
 Resids_msis86) = Read_GEODYN_func(arc1,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_msis86,
                            True,
                            2003,
                            'SLR')

(AdjustedParams_dtm87,
 Trajectory_dtm87,
 Density_dtm87,
 Resid_dtm87s) = Read_GEODYN_func(arc1,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_dtm87,
                            True,
                            2003,
                            'SLR')

(AdjustedParams_jaachia,
 Trajectory_jaachia,
 Density_jaachia,
 Resids_jaachia) = Read_GEODYN_func(arc1,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_jaachia71,
                            True,
                            2003,
                            'SLR')
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

There are some weird things here for Gen.Acc. Not done yet.
There are some weird things here for Gen.Acc. Not done yet.
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/dtm87/dtm87_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/dtm87/dtm87_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/dtm87/dtm87_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

There are some weird things here for Gen.Acc. Not done yet.
There are some weird things here for Gen.Acc. Not done yet.
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/jaachia71/jaachia71_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/jaachia71/jaachia71_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/jaachia71/jaachia71_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

There are some weird things here for Gen.Acc. Not done yet.
There are some weird things here for Gen.Acc. Not done yet.
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

Drag Coefficient Comparison

[4]:
import plotly.graph_objects as go
import numpy as np

labels = list(AdjustedParams_msis86[5][SAT_ID]['0CD'].keys())

Cd_plot = []
for i in AdjustedParams_msis86.keys():
    Cd_plot.append(float(AdjustedParams_msis86[i][SAT_ID]['0CD'][labels[0]]['CURRENT_VALUE']))
fig = go.Figure(data=[go.Scatter(x=list(AdjustedParams_msis86.keys()),
                                   y=Cd_plot,
                                 name = 'MSIS 86',
                               mode='markers',
                                    marker=dict(
                                    size=10,
                                    ),
                                   ),
                                   ],
                                   )

Cd_plot = []
for i in AdjustedParams_dtm87.keys():
    Cd_plot.append(float(AdjustedParams_dtm87[i][SAT_ID]['0CD'][labels[0]]['CURRENT_VALUE']))
fig.add_trace(go.Scatter(x=list(AdjustedParams_dtm87.keys()),
                                   y=Cd_plot,
                                      name = 'DTM 87',
                                   mode='markers',
                                  marker=dict(
                                    size=10,
                                    ),
                                   ),
                                 )
Cd_plot = []
for i in AdjustedParams_jaachia.keys():
    Cd_plot.append(float(AdjustedParams_jaachia[i][SAT_ID]['0CD'][labels[0]]['CURRENT_VALUE']))
fig.add_trace(go.Scatter(x=list(AdjustedParams_jaachia.keys()),
                                   y=Cd_plot,
                                 name = 'Jaachia 71',
                               mode='markers',
                                    marker=dict(
                                    size=10,
                                    ),
                                   ),
                                 )

fig.update_layout(
    title="Drag Coefficient vs Iteration",
    yaxis_title=r"$C_d \,\,\text{(Unitless)}$",
    xaxis_title="Iteration Number",
    )

fig.show()


The drag coefficient is time dependent

[5]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams_msis86[5][SAT_ID]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=("Time Dependent Drag Coefficient (last iter)",
                    ))
val_list = []
for i in AdjustedParams_msis86[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_msis86[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'MSIS 86',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_dtm87[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_dtm87[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_jaachia[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_jaachia[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'Jaachia 71',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Date", row=1, col=1)

iplot(fig)

Compare density along the orbit

[6]:
data_nums = 1000

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][:data_nums],
                                 y=Density_msis86['rho (kg/m**3)'][:data_nums],
                                    name='MSIS 86',
                                     mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=Density_dtm87['Date'][:data_nums],
                                 y=Density_dtm87['rho (kg/m**3)'][:data_nums],
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=2,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][:data_nums],
                         y=Density_jaachia['rho (kg/m**3)'][:data_nums],
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )

fig.update_layout(
    title="Density along Starlette Orbit (for 4 minutes of first day)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)








[7]:
data_nums = 100
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][::data_nums],
                                 y=Density_msis86['rho (kg/m**3)'][::data_nums],
                                 name = 'MSIS 86',
                                 mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.add_trace(go.Scatter(x=Density_dtm87['Date'][::data_nums],
                                 y=Density_dtm87['rho (kg/m**3)'][::data_nums],
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=2,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][::data_nums],
                         y=Density_jaachia['rho (kg/m**3)'][::data_nums],
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )
fig.update_layout(
    title="Density along Starlette Orbit (Every 100th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)






Compare residuals

[8]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Resids_msis86['Date'],
                                 y=Resids_msis86['Residual'].values.astype(float),
                                 name ='MSIS 86' ,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=Resid_dtm87s['Date'],
                                 y=Resid_dtm87s['Residual'].values.astype(float),
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=3,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=Resids_jaachia['Date'],
                         y=Resids_jaachia['Residual'].values.astype(float),
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )

iplot(fig)





Compare Trajectories:

[9]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Trajectory_msis86['Date'],
                                 y=Trajectory_msis86['Z'].values.astype(float),
                                 name ='MSIS 86' ,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=Trajectory_dtm87['Date'],
                         y=Trajectory_dtm87['Z'].values.astype(float),
                         name= 'DTM 87',
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )
fig.add_trace(go.Scatter(x=Trajectory_jaachia['Date'],
                         y=Trajectory_jaachia['Z'].values.astype(float),
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )

iplot(fig)





[ ]: